#ifndef RTNORM_HPP
#define RTNORM_HPP

#include "../Random.hpp"

//Implementation of Robert's truncated Gaussian (1995)

#define pINF 10000

class rtnorm
{
	public:
	rtnorm()
	{
		boost::random::uniform_01<> U;
		runif= new RandomG::Random<boost::random::uniform_01<> >(U);
	}
	~rtnorm()
	{
		delete runif;
	}
	double operator()(double m,double mm,double mp,double s)
	{
		double res=0;
		if(-mm==pINF)
		{
						
			mp=(mp-m)/s;			
			mm=-mp;
			double alpha=0.5*(mm+sqrt(pow(mm,2)+4));
			int i=0;
			while(i==0)
			{
				double u=(*runif)();
				double x=(-1/alpha)*log(1-u)+mm;
				double rho=exp(-0.5*pow((x-alpha),2));
				double u2=(*runif)();
				if(u2<=rho)
				{
					res=(-x*s+m);
					i=1;
				}
			}
			
		}else if(mp==pINF){
			mm=(mm-m)/s;
			double alpha=0.5*(mm+sqrt(pow(mm,2)+4));
			int i=0;
			while(i==0)
			{
				double u=(*runif)();
				double x=(-1/alpha)*log(1-u)+mm;
				double rho=exp(-0.5*pow((x-alpha),2));
				double u2=(*runif)();
				if(u2<=rho)
				{
					res=x*s+m;
					i=1;
				}
			}
		}else{
			double i=0;
			mp=mp-m;
			mm=mm-m;
			while(i==0)
			{
				double z=(mp-mm)*(*runif)()+mm;
				double r=rho(z,mm,mp);	
				double u=(*runif)();
				if(u<=r)
				{
					i=1;
					res=z+m;				
				}
			}
		}
		return res;
	}
	double rho(double z, double mm, double mp)
	{
		double res=0;
		if(0>mp)
		{
			res=exp((pow(mp,2)-pow(z,2))/2);
		}else if(mm>0){
			res=exp((pow(mm,2)-pow(z,2))/2);
		}else{
			res=exp(-(double)pow(z,2)/2);
		}
		return res;
	}

	private:
	RandomG::Random<boost::random::uniform_01<> > *runif;
};

#endif

